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A  feature  of  many  ill-posed  inverse  problems  is  that  the  Hessian  operator  of  the  data 
misfit  functional  is  a  compact  operator  with  rapidly  decaying  eigenvalues.  This  is  a 
manifestation  of  the  typically  sparse  observations,  which  are  informative  about  a  limited 
number  of  modes  of  the  infinite  dimensional  held  we  seek  to  infer.  The  Hessian  operator 
(and  its  finite  dimensional  discretization)  play  an  important  role  in  the  analysis  and 
solution  of  the  inverse  problem.  In  particular,  the  spectrum  of  the  Hessian  at  the 
solution  of  the  inverse  problem  determines  the  degree  of  ill-posedness  and  provides 
intuition  on  the  construction  of  appropriate  regularization  strategies.  This  has  been 
observed,  analyzed,  and  exploited  in  several  applications  including  shape  optimization 
[1,  2]  and  inverse  wave  propagation  [3,  4,  5],  to  name  a  few. 

Moreover,  solution  of  the  inverse  problem  by  the  gold  standard  iterative  method — 
Newton’s  method — requires  “inversion”  of  the  Hessian  at  each  iteration.  Compactness 
of  the  Hessian  of  the  data  misfit  functional  accompanied  by  sufficiently  fast  eigenvalue 
decay  permits  a  low  rank  approximation,  which  in  turn  facilitates  rapid  inversion  or 
preconditioning  of  the  regularized  Hessian  [3,  6].  Alternatively,  solution  of  the  linear 
system  arising  at  each  Newton  iteration  by  a  conjugate  gradient  method  can  be  very  fast 
if  the  data  misfit  Hessian  is  compact  with  rapidly  decaying  eigenvalues  and  the  conjugate 
gradient  iteration  is  preconditioned  by  the  regularization  operator  [7].  Finally,  under  a 
Gaussian  approximation  to  the  Bayesian  solution  of  the  inverse  problem,  the  covariance 
of  the  posterior  probability  distribution  is  given  by  the  inverse  of  the  Hessian  of  the 
negative  log  likelihood  function.  For  Gaussian  data  noise  and  model  error,  this  Hessian 
is  given  by  an  appropriately  weighted  Hessian  of  the  data  misfit  operator,  e.g.,  [8].  Here 
again,  exploiting  the  low-rank  character  of  the  data  misfit  component  of  the  Hessian  is 
critical  for  rapidly  approximating  its  inverse,  and  hence  the  uncertainty  in  the  inverse 
solution  [4,  5,  9,  10]. 

In  all  of  the  cases  described  above,  compactness  of  the  data  misfit  Hessian  is  a 
critical  feature  that  enables  fast  solution  of  the  inverse  problem,  scalability  of  solvers  to 
high  dimensions,  and  estimation  of  uncertainty  in  the  solution.  With  this  motivation, 
here  we  analyze  the  Hessian  operator  for  inverse  medium  acoustic  scattering  problems, 
and  study  its  compactness.  Our  analysis  is  based  on  an  integral  equation  formulation 
of  the  Helmholtz  equation,  adjoint  methods,  and  compact  embeddings  in  Holder  and 
Sobolev  spaces.  These  tools  allow  us  to  analyze  the  shape  Hessian  in  detail. 

The  remainder  of  the  paper  is  organized  as  follows.  Section  2  briefly  derives 
and  formulates  forward  and  inverse  acoustic  scattering  problems  due  to  bounded 
inhomogeneity.  We  then  derive  the  Hessian  for  the  inverse  problem  in  Section  3.  Section 
4  justifies  the  Hessian  derivation  by  studying  the  well-posedness  of  the  (incremental) 
forward  and  (incremental)  adjoint  equations,  and  the  regularity  of  their  solutions.  Next, 
we  analyze  the  Hessian  in  Holder  spaces  in  Section  5,  and  then  extend  the  analysis  to 
Sobolev  spaces  in  Section  6.  In  order  to  validate  our  theoretical  developments,  we 
provide  numerical  examples  in  Section  7.  Finally,  the  conclusions  of  the  paper  are 
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2.  Forward  and  Inverse  Medium  Problems  For  Acoustic  Scattering 


In  this  section,  we  briefly  discuss  forward  acoustic  scattering  problems  due  to  bounded 
inhomogeneity  and  the  corresponding  inverse  problems.  Since  both  forward  and  inverse 
medium  problems  can  be  found  elsewhere  [11],  our  attention  is  to  introduce  necessary 
notations  that  will  be  used  in  our  later  derivations  and  analysis  starting  from  Section  3. 

The  scattering  of  time  harmonic  incident  wave  due  to  bounded  inhomogeneity  can 
be  shown  to  be  governed  by  the  following  Helmholtz  equation  [11]: 

V2!/  +  k2nU  =  (1  -  n)k2Uic,  in  Md,  (la) 


hm  r^2  I  ^ 
or 


ikU  )  =  0, 


r  = 


(lb) 


where  Ulc  is  the  incident  wave  that  satisfies  the  Helmholtz  equation  V2f/*c  +  kUlc  =  0, 
k  >  0  the  wave  number,  n  >  0  the  refractive  index  which  is  assumed  to  be  1  for  the  free 
space,  d  G  {2,3}  the  dimension  of  the  background  space,  i2  =  —1,  and  U  the  scattered 
field.  The  radiation  condition  (lb)  is  assumed  to  be  valid  uniformly  in  all  directions 
jj^jl  with  x  denoting  the  vector  of  spatial  coordinates.  To  the  rest  of  the  paper,  the 
inhomogeneity  is  assumed  to  be  bounded,  that  is,  there  exists  some  sufficiently  large 
a  >  0  such  that  n(x)  =  1,  V  ||x||  >  a.  In  other  words,  q  =  1  —  n  has  compact  support  in 
Md. 


For  the  forward  problem,  n  is  given  and  we  solve  the  forward  equations  (la)-(lb) 
for  the  scattered  field  U.  For  the  inverse  problem,  on  the  other  hand,  given  observation 
data  Uobs  over  some  compact  subset  klobs  C  Md,  we  are  asked  to  infer  the  distribution 
of  the  refractive  index  n.  One  way  to  solve  the  inverse  problem  is  to  cast  it  into  the 
following  PDE-constrained  optimization  problem: 


min  J  =  [  K(x)  I U  -  Uobsf  dQ,  (2) 

q  JRd 

subject  to  the  forward  equations  (la)-(lb).  Here,  A'(x)  denotes  the  observation  operator 
whose  support  is  Qobs.  In  order  to  cover  several  interesting  observation  operators,  Qobs 
is  allowed  to  be  quite  general  in  this  paper.  In  particular,  it  could  be  a  closed  subset  in 
or  a  relative  closed  subset  of  a  manifold  in  Md.  For  example  in  M3,  Q°bs  could  be  a 
closed  arc,  or  a  closed  curved,  or  a  closed  subset  of  a  two  dimensional  manifold,  or  some 
two  dimensional  manifold.  For  convenience,  we  identify 

IC(p  =  I  Kip  dl  l  —  j  (p(y)dy. 
jRd  Jn°bs 

h  (  hi  Noba 

We  also  permit  pointwise  observation  in  our  analysis,  i.e.,  Slobs  =  {xj  }^._1  ,  and  in  this 
case  we  identify 
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In  this  section,  we  derive  the  gradient  and  Hessian  using  a  reduced  space  approach, 
and  the  justification  for  our  derivations  is  provided  in  Section  4.  We  begin  with  a 
useful  observation  on  the  radiation  condition.  Since  the  radiation  condition  (lb)  is  valid 
uniformly  in  all  directions  p^,  we  rewrite  the  radiation  condition  as 

^ -  —  ikU  =  (p(r)  =  o  (r(1-d)/2)  , 
or 

where  r  is  the  radius  of  a  sufficiently  large  circle  Too. 

It  can  be  seen  that  the  cost  functional  (2)  is  real- valued  while  the  constraints  (la)- 
(lb)  are  complex-valued.  Consequently,  the  usual  Lagrangian  approach  will  not  make 
sense  and  care  must  be  taken.  Following  Kreutz-Delgado  [12],  we  define  the  Lagrangian 
as 


C  =  J+[  u[V2U  +  k2nU 

JRd 


V2U  +  k2nU-k2(  1 


k2(l  -n)Uic]  dQ  + 


ds , 


ds 


where  the  overline,  when  acting  on  forward  and  adjoint  states  (and  their  variations), 
denotes  the  complex  conjugate. 

Taking  the  first  variation  of  the  Lagrangian  with  respect  to  u,  ur  in  the  directions 
u,  ur  and  arguing  that  the  variations  u,  ur  are  arbitrary  yield  the  forward  equations 
(la)-  (lb). 

Now  taking  the  first  variation  of  the  Lagrangian  with  respect  to  U  in  the  direction 
U  and  arguing  that  the  variation  U  is  arbitrary  yield  the  following  adjoint  equations: 


V2w  +  k2nu  =  -K  ( U  -  Uobs )  , 

lim  r^_1)/2  ( ^  +  ikv\  =  0, 

? — »oo  V  or  ) 


m 


r  =  x 


(4a) 

(4b) 


and 


ur  =  —u  on  Too. 

If  we  eliminate  ur:  the  Lagrangian  now  becomes 

C  =  J+[  u  [V2U  +  k2nU  -  k2{  1  -  n)Uic]  dSl  - 


u 


dU 

dr 


—  ikU  —  ip)  ds 


u 


V2U  +  k2nU-  k2(l-n)Uic  dU-J  u  -  ikU  -<p  )  ds.  (5) 


The  gradient  of  the  cost  function  acting  in  the  direction  n  is  simply  the  variation  of  the 
Lagrangian  with  respect  to  q  in  the  direction  h,  i.e., 


DJ(q;h)  =  -k2  u(U  +  Uic)  +u[U  +  U 


n  dM. 


(6) 
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For  the  sake  of  convenience  in  deriving  the  Hessian,  the  forward  and  adjoint 
equations  are  best  expressed  in  the  weak  form.  As  a  direct  consequence  of  the  above 
variational  calculus  steps,  the  forward  equation  in  the  weak  form  reads 


S(q,U)=  u  [V2C/  +  k2nU  -  k2(l  -  n)Uic]  dtt 


— - ikU  —  (p  1  ds  =  0,  Vu. 

or 


Similarly,  the  adjoint  equation  in  the  weak  form  is  given  by 


(7) 


A  (q,  U,  u ) 


u  +  k2nu 


+  K  (U  —  Uobs)]  dn 


iku 


ds 


0,  vu. 


(8) 


Next,  the  (reduced)  Hessian  acting  in  the  directions  n  and  n  is  obtained  by  simply 
taking  the  first  variation  of  the  gradient  DJ  (q,  h)  with  respect  to  q,  U  and  u  in  the 
directions  n,  U  and  u,  i.e., 


D2J(q;h,n) 


u(U  +  Uic )  +u(u  +  Ulc) 


+  uU  +  uU 


ndkl. 


(9) 


As  mentioned  at  the  beginning  of  this  section,  the  reduced  space  approach  is  employed, 
and  hence  the  variations  U  and  u  can  not  be  arbitrary.  In  fact,  they  are  only  admissible 
if  the  forward  and  adjoint  equations  are  satisfied.  As  a  direct  consequence,  the  first 
variations  of  S  (q,  U )  and  A  (q,  U,  u )  must  vanish,  that  is,  U  is  the  solution  of  the 
following  incremental  forward  equation: 


V2U  +  k2nU  -  k2n  (U  +  Uic ) 


dfl 


ds 


0,  Vu, 


and  e  is  the  solution  of  the  following  incremental  adjoint  equation: 


(10) 


u 


+  k2nu  —  k2hu  +  KU 


dM 


!  U  t^  +  iku)  ds  =  0,  VC/. 

'Too  \ 


(11) 


Consequently,  the  corresponding  strong  form  of  the  incremental  forward  equation  is 


V2CC  +  k2nU  =  k2n  ( U  +  Uic )  , 


lirn  r(d"1)/2 

r— >oo 


0, 


in  Rd,  (12a) 

r  =  ||x|| ,  (12b) 
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in  Md,  (13a) 

r  =  ||x||  .  (13b) 

Next,  we  need  to  convert  the  Hessian  in  (9)  into  a  form  that  is  convenient  for  our  later 

analysis.  The  first  step  is  to  replace  h  by  n  and  choose  u  =  u{n )  in  the  incremental 

forward  equation  (10).  In  the  second  step,  we  take  U  =  U{n)  in  the  incremental  adjoint 
equation  (11).  The  last  step  is  to  subtract  the  resulting  incremental  forward  equation 
from  the  complex  conjugate  of  the  resulting  incremental  adjoint  equation.  After  some 
simple  integration  by  parts  and  cancellations,  we  obtain 


V2w  +  k2nu  =  k2nu  —  KU, 


du 


lim  r (d  1)/2  (  — — b  iku  )  =  0, 
or 


/  k2nu(n)(U  +  Uic)dn=  k2unU(n)dn~  KU(h)U(h)dn.  (14) 

jRd  JRd  JRd 

Combining  (9)  and  (14)  gives  the  desired  form  of  the  Hessian  as 


D2  J (g;  h,  n)  =  /  K  U (n)U (h)  +  U (h)U (h) 


dM 


uU  (n)  +  uU  (n) 


T-Li{q\h,n) 

n  dll  +  k 2 


uU  (h)  +  uU  (n) 


n  dll 


(15) 


4.  Regularity  of  the  forward  and  adjoint  solutions 


In  this  section  we  are  going  to  justify  what  we  have  done  in  Section  3  by  studying  the 
wcll-posedness  of  the  forward  and  adjoint  equations,  and  the  regularity  of  their  solutions. 
For  sufficiently  smooth  inhomogeneity,  the  solutions  turn  out  to  be  classical  by  using 
an  integral  equation  method,  as  we  shall  show. 

First  we  introduce  the  following  standard  volume  potentials  (also  known  as  Newton 
potentials)  [13,  11]: 

w(x)  =  7V(x)  =  f  *h(x,  y)</?(y)  dy,  x  G  (16) 

J  RN 


where  is  the  fundamental  solution  of  the  (incremental)  forward  equation(s)  defined  as 


$(x,y) 


l^o1  (x  -  y) 

gifellx  — y|| 

4vr||x-y|| 


N 

N 


2 

3  ’ 


or  the  fundamental  solution  of  the  (incremental)  adjoint  solution(s): 


-IH2(, 

e-»fc||x-y|| 
4tt||x — y  || 


y)  n 

N 


2 

3 


$(x,y) 
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Next  we  denote  Cm,a  (Rd)  the  space  of  m- times  differentiable  functions  whose  mth 
derivative  is  Holder  continuous  with  exponent  a  in  Rd,  and  C™’a  (Rd)  a  subspace 
of  Cm,a  (Rd)  consisting  of  functions  with  compact  support.  The  following  mapping 
properties  of  the  volume  potential  are  important  in  the  sequel. 

Lemma  1  Let  ip  G  C0  (Md)  D  where  a  G  (0,1]  and  me  NU  {0}.  Then 

w  e  Cm+2’a(fl),  and  ||w||Cm+2,a(Q)  <  c  where  supp(ip)  C  f 1C  Rd. 

Furthermore,  it  holds  that  T  is  compact  in  Cp,^(fl)  for  m  +  a  <  p  +  (3  <m  +  2  +  a. 

Proof.  We  proceed  by  induction.  The  cases  m  =  0, 1  are  proved  and  discussed  in 
[13,  11].  Since  differentiation  and  integration  can  be  interchanged  [13,  11],  the  partial 
derivative  in  Xj  direction  reads 


DXjw{x)  =  /  <p(y)Dx$  (x,  y)  dfl  =  -  (p(y)Dy$  (x,y)  dy 


=  -  $  (x,  y)  (f(y)rij  ds+  (x,  y)  Dy.(p( y)  dy, 

Jon  Jn 


(17) 


where  nj  denotes  the  jth  component  of  the  normal  vector  of  dfl,  and  we  have  used  the 
fact  that  DXj  &  (x,  y)  =  —  Dyj <f>  (x,  y)  and  (p(y)  =  0  on  dll.  Since  Dy.tp( y)  e  C0  (Md)  D 
Cm-i,a(Md),  we  conclude,  by  the  induction  hypothesis,  that  DXjw(x )  G  Cm+1,a  (ft). 
This  implies  w(x)  G  Cm+2,a  (ft). 

The  second  assertion  is  readily  proved  by  the  induction  hypothesis 


and  the  definition  of  Holder  norms  [14,  11],  The  third  assertion  is  trivial  due  to  the 
compact  embeddings  in  Holder  spaces  [14].  □ 

In  order  to  use  Lemma  1  and  the  Riesz-Fredholm  theory  [15,  16]  to  study  the 
well-posedness  of  the  forward  and  adjoint  equations,  we  first  recall  the  following  Green 
formula  [11]  for  u  G  C2(fl)  D  C(fl) 


ux  = 


It 90  L 


du  <9$ 
“  Udy_ 


ds  -  I  4>  (V2m  +  k2u)  dfl, 


where  n  denotes  the  unit  outward  normal  vector  of  dfl.  Denote 


T[q]p(x)  =  k2  [  <&(x,  y)q(y)p(y)  dy,  xGO, 

Jn 

and  /  as  the  identity  operator.  We  now  have  the  following  integral  representations  for 
the  forward  and  the  adjoint  equations. 

Theorem  1  Let  ft  C  Rd,  mid  q,h  G  Gon’“(fl),  where  a  G  (0, 1]  and  m  G  MU  {0}.  In 
addition,  let  p  and  /3  satisfy  m  +  a  <  p  +  (3  <  m  +  2  +  a.  The  forward,  incremental 
forward,  adjoint,  and  incremental  adjoint  equations  are  well-posed  in  the  sense  that  they 
are  equivalent  to  the  following  Lippmann- Schwinger- type  integral  equations: 
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i)  The  forward  integral  equation 


{I  +  T[q\)  t/(x)  =  -T[q]Uicfx)  (18) 

has  a  unique  solution  U  in  Cp,/3(£2)  and  the  solution  depends  continuously  on  the 
data. 

ii)  The  incremental  forward  integral  equation 

(/  +  T[q\)  *7(x)  =  T[h]  (U  +  Uic )  (x)  (19) 

has  a  unique  solution  U  in  Cp,/3(£2)  and  the  solution  depends  continuously  on  the 
data. 

in)  The  adjoint  integral  equation 

(/  +  T[9])u(x)=  [  Q(x,y)K  (U  —  Uobs)  (y)dy  (20) 

Jn°bs 

has  a  unique  solution  in  Cm,a  (£1)  and  the  solution  depends  continuously  on  the 
data. 

iv)  The  incremental  adjoint  integral  equation 

(/  +  T[q\)  «(x)  =  T[n]u(x)  +  [  <F(x,  y )KU (y)  dy  (21) 

Jn°bs 

has  a  unique  solution  in  Cm,a  (£2)  and  the  solution  depends  continuously  on  the 
data. 

Proof.  The  equivalence  for  the  forward  equation  can  be  found  in  [11,  17]  for  the 
case  q  e  Cq(M3).  The  generalization  to  the  case  q  G  is  straightforward.  By  the 

same  token,  we  have  the  equivalence  for  other  equations.  This  suggests  that  we  need  to 
study  only  the  well-posedness  of  the  integral  equations. 

By  Lemma  1,  T  is  compact  in  Cp,/3(£2),  and  T[g]7Rc(x)  e  Cp,/3(f2)  due  to 
the  analyticity  of  the  incident  wave  Ulc.  The  Riesz-Fredholm  theory  [15]  therefore 
applies,  and  the  forward  integral  equation  has  a  unique  solution  in  Cp,^(Tl).  Moreover, 
(/  +  T[q\Yl  is  bounded,  that  is,  the  solution  depends  continuously  on  T  [g]  Ulc (x)  in 
the  Cp,P- norm.  The  proof  for  incremental  forward  solution  (19)  follows  the  same  line 
by  observing  that  the  right  side  of  (19)  belongs  to  Cp,/3(£2)  from  Lemma  1. 

As  for  the  adjoint  equation,  owing  to  Plobs  fl  supp  (q)  =  0  and  the  analyticity  of 
$(x,  y),  the  right  side  of  (20)  is  certainly  a  function  in  Cp,l3(Q)  (in  fact  it  is  analytic  in 
Md\supp  (£2o6s)),  and  again  the  Riesz-Fredholm  theory  gives  the  desired  results.  Finally, 
for  the  incremental  adjoint  integral  equation,  observe  that  the  first  term  on  the  right 
side  of  (21)  belongs  to  Cp,^{Pl)  and  the  second  term  is  analytic  on  £2.  As  a  result,  the 
right  hand  side  of  equation  (21)  is  function  in  CP,^(Q).  The  conclusions  are  now  readily 
verified  by  the  Riesz-Fredholm  theory.  □ 

We  are  now  in  the  position  to  justify  our  derivations  of  gradient  and  Hessian  in 
Section  3. 
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Theorem  2  Let  £7  C  Rd  be  a  bounded  domain.  Assume  that  q,  h  and  h  belong  to 
Cq  (£7).  Then,  the  cost  functional  (2)  is  twice  continuously  Frechet  differentiable;  and 
hence  the  gradient  (6)  and  Hessian  (15)  are  well-defined. 

Proof.  First,  observe  that  we  have  used  Gateaux  derivatives  to  derive  the  gradient 
and  Hessian  in  Section  3.  Now  it  is  evident  that  both  DJ (q\  h)  and  D'2J (q;  h,  h)  are 
linear  and  continuous  with  respect  to  h  (and  ft)  since  U,  u,  U,  and  u  belong  to  C2,a  (£7) 
by  Theorem  1.  Moreover,  continuous  dependence  on  q  of  U  from  Theorem  1  implies 
the  continuous  dependence  on  q  of  u,  U,  and  u,  which  in  turn  implies  the  continuity  of 
DJ  (q;  h)  and  D2J  (q;  h,  h)  with  respect  to  q.  Hence,  a  classical  result  on  sufficiency  for 
Frechet  derivative  [18]  ends  the  proof.  □ 


5.  Analysis  of  the  Hessian  in  Holder  spaces 

In  this  section  we  study  the  behavior  of  the  Hessian  at  a  fixed  refractive  distribution 
n,  i.e.,  q  —  1  —  n  G  Cff’a  (£7).  Unlike  the  shape  Hessian  which  is  only  compact  at  the 
optimal  solution  as  we  have  analyzed  in  the  first  part  of  this  work  [19],  the  Hessian 
of  the  inverse  medium  scattering  problem  turns  out  to  be  compact  for  all  q  as  we 
shall  show.  For  concreteness,  we  restricted  ourselves  to  two  exemplary  cases  of  the 
observation  operator,  namely,  the  observation  is  everywhere  on  a  compact  subset  Llobs 
having  non-trivial  r-dimensional  Lesbegue  measure  for  some  1  <  r  <  d  (we  call  this 
case  as  continuous  observation)  and  pointwise  observation  £7obs  =  jx^  ]■  . 

From  Theorem  1,  observe  that  the  incremental  forward  solution  U  can  be  identified 
as  the  following  operator  composition: 

U  :  C™’°  (£7)  3h^  U(n)  =  (/  +  T[g])_1  T[n]  ( U  +  Uic)  e  C^(i 7), 

which  is  compact  since  it  is  the  composition  of  the  continuous  operator  (/  +  T[g])  1 
(owing  to  the  Riesz-Fredholm  theory)  and  the  compact  operator  T[n]  ( U  +  Ulc )  (due  to 

Lemma  1).  As  a  result,  U (h)  is  still  a  compact  operator  since  the  restricting  to  Dobs 

a.obs 

is  a  continuous  operation. 

If  the  observation  is  continuous,  the  Gauss-Newton  part  of  the  Hessian,  namely 
Hi  (n;n,  n),  can  be  now  rewritten  as 


Hi(q;  h,  h)  =  2^  (  U (h),  U (h 


L2(flobs) 


=  2S?  [U*U(h),n 


L2(fl°bs) 


where  the  real  operator  &  extracts  the  real  part  of  its  argument,  and  (•)*  denotes  the 
adjoint  operator.  In  this  form,  H\(n)  is  evidently  compact  due  to  the  compactness  of 
U  (n) 

noba 

If,  on  the  other  hand,  the  observation  is  pointwise,  then  the  evaluation  of  U (h)  at 
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x.°jhs  can  be  written  as 

U(n)  (xf )  =  {-k\^,U(h))m^  +  + 

=  (k%  (u  +  u*)  -  in  (fc2^)  ,n)  LHn) 

=  (H/j,  h)L 2^Q)  , 

where  =  $(x°6s,y)  and  =  k2&j  (U  +  Ulc )  —  U*  (k2q$j).  In  this  case,  the  Gauss- 
Newton  part  Hi(n-,h,  h)  reads 


( . 

( Nobs  \ 

\ 

Hi(q-,  h,  h)  =  2& 

E^VIv^ 

1  ,fl 

V 

V  ) 

L2m  J 

which  shows  that  the  dimension  of  the  range  of  H\ (q)  is  at  most  Nobs.  Consequently, 
Hi  (q)  is  a  compact  operator.  We  summarize  the  above  result  on  the  compactness 
of  Hi{q)  in  the  following  theorem  which  is  valid  for  both  continuous  and  pointwise 
observation  cases. 

Theorem  3  Hi{q),  as  a  continuous  bilinear  form  on  C™’a  (Rd)  x  C™'a(Rd),  is  a 
compact  operator. 

The  analysis  of  H-iiq)  is  somewhat  easier  as  we  shall  now  show. 

Theorem  4  H2(q),  as  a  continuous  bilinear  form  on  Cfl'°  x  Cfl'°  (Md),  is  a 
compact  operator. 

Proof.  Rewrite  H2(q',h,h)  as 


H2(q ; 


=  2  k2& 


uU (h)n  +  uU (n)n  dfl  =  2 k2S?  ( U*  ( un )  +  uU (h),  n 


L2(0) 


We  conclude  that  H2  ( q )  is  compact  by  the  following  three  observations.  First,  the 
incremental  forward  solution  U  can  be  identified  as  a  compact  operator  in  C™’0  (fl)  as 
discussed  above.  Second,  multiplication  by  u  G  Gq1’0  (G)  is  a  continuous  operation  (see 
[19]  for  example).  Third,  the  sum  of  two  compact  operators  is  again  compact.  □ 

We  close  this  section  by  observing  that  the  full  Hessian  is  the  difference  of  two 
compact  operators,  it  is  therefore  compact  as  well. 


6.  Analysis  of  the  Hessian  in  Sobolev  spaces 

Similar  to  the  first  part  of  our  work  [19],  we  shall  extend  the  analysis  in  Holder  spaces 
to  Sobolev  spaces.  A  result  similar  to  that  of  Lemma  1  is  now  stated. 

Lemma  2  Assume  that  Lp  is  bounded  and  integrable,  Q  C  Rd  is  a  bounded  domain, 
and  supp(tp)  C  fb  Then  T  defined  in  (16)  maps  Hrn  (G)  continuously  to  Hm+ 2  (f!)  for 
rn  G  N  U  {0}. 
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Proof.  We  proceed  by  induction.  Case  rn  —  0  has  been  already  proved  in  [11], 
Now  assume  that  the  assertion  holds  for  m  —  1,  and  we  need  to  show  that  it  also  holds 
for  m.  Since  boundedness  and  integrability  of  ip  enable  integration  and  differentiation 
interchange  [20],  (17)  holds.  By  the  induction  hypothesis  DXjw(x.)  G  H'n+1  (Q),  and 
this  implies  w  (x)  =  Tip  (x)  G  Hm+2  (fl).  □ 

By  compact  embeddings  in  Sobolev  spaces  [21,  22],  one  can  see  that  T  is  compact 
in  Hs  (fl)  for  m  <  s  <  m  +  2.  This  fact  is  used  to  prove  the  following  compactness  of 
the  Hessian  in  Sobolev  spaces. 

Theorem  5  Let  q  be  bounded,  integrable,  and  q  G  Hm  ( O) ,  where  f 1  is  a  bounded 
domain,  and  supp(q)  C  fl  C  Rd.  Then,  the  Hessian,  H{q)  =  Hi(q)  —  Hi^q),  is  a 
compact  operator  in  Hm  (fl) . 

Proof.  The  proof  follows  the  same  line  as  in  Section  5  by  using  Lemma  2,  and  hence 
omitted.  □ 

7.  Numerical  results 

In  this  section,  we  numerically  compute  the  eigenvalues  of  the  shape  Hessian  (9) 
to  validate  our  theoretical  developments  in  Sections  5  and  6.  For  the  purpose  of 
demonstration,  it  is  sufficient  to  consider  two-dimensional  problems  for  which  we  can 
use  an  efficient  coupled  finite  element  and  boundary  integral  equation  approach.  The 
detailed  description  of  our  coupling  strategy  is  now  presented. 

7.1.  Forward  scattering  problem 

We  decompose  the  forward  problem  into  two  sub-problems,  namely,  the  interior  sub¬ 
problem  given  by 

V2f/in  +  k2nUin  =  (1  -  n)k2Uic , 
dUin 

—  +ikUin  =  if, 

on 

and  the  exterior  sub-problem  given  by 


in  Rd  \  hi,  (23a) 

on  d H,  (23b) 

(23c) 

where  if  is  the  unknown  coupling  function.  It  is  easy  to  show  that  the  weak  formulation 
of  the  interior  problem  (22)  reads 

f  Vf/m  -Vw  dPl  —  k2  f  nUmv  dtt+ik  f  Umv  ds  =  ik  f  ifvds  —  k2  f  qUlcv  d,Q,  (24) 
Jn  Jn  Jon  Jon  Jn 


in  hi, 

(22a) 

on  <9f 1, 

(22b) 
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and  that  the  boundary  integral  equation  formulation  of  the  exterior  problem  (23)  reads 


(. I-D -  ikS )  Uex  =  -Sit), 


(25) 


with  the  representation 


t/ex(x)  =  /  Uex(y) 
■Jon 


d$  (x,  y) 

dn 


+  ik&  (x,  y) 


ds( y)  -  /  $(x,y)^(y)ds(y), 


*  on 


and  with  S,D  as  the  following  standard  surface  single  and  double  layer  potentials  [11]: 


SV(x)  =  2/  $(x,  y)<p(y)  ds(y), 

J  80. 

D<p{x)  =  2  f  ~ ds(y), 
J  do  9n(y ) 


x  G  <90, 
x  G  (90. 


The  interior  and  exterior  solutions  are  matched  by  satisfying  the  following 
continuity  condition  at  the  interface  <90: 


U in  =  Uex ,  on  on, 


(26) 


which,  together  with  (22)  and  (23),  implies  the  continuity  in  the  normal  derivative,  i.e., 


dUin 

dn 


8Uex 

dn 


on  (90. 


Inspired  by  the  coupling  approach  in  [23],  we  choose  to  use  the  finite  element  method 
(FEM)  for  solving  the  interior  problem  (24),  while  we  solve  the  exterior  boundary 
integral  equation  (25)  using  a  Nystrom  method  [11,  19].  Now,  the  nature  of  the  coupling 
is  implicit,  that  is,  in  order  to  solve  for  Um  and  Uex,  the  availability  of  ip  is  required.  On 
the  other  hand,  in  order  to  solve  (26)  for  ip,  one  has  to  supply  U™  and  Uex.  Moreover, 
matching  the  finite  element  and  Nystrom  methods  may  not  be  trivial  since  the  finite 
element  solution  is  defined  variationally,  while  the  Nystrom  solution  is  pointwise  in 
nature.  We  adopt  a  simple  decoupling  approach  due  to  Kirsch  and  Monk  [23]  in  which 
ip  is  represented  by  trigonometric  polynomials  of  order  M.  i.e., 


M—l 

ip=  ^2  aj&j,  (Pj  =  ezjt,  te[ 0,2vr], 

j=—M 


With  this  representation,  one  can  solve  the  interior  and  exterior  problems  independently 
for  each  basis  function  (fij.  Then,  the  unknown  coefficients  ay  can  be  solved  for  by 
employing  a  Galerkin  projection  on  (26),  i.e., 

[  (; Uin  -  Uex)  4>j  ds  =  0,  j  =  -M, . . . ,  M  -  1, 

JdO 


where  dj  denotes  the  complex  conjugate  of  dj- 
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In  this  paper,  we  use  the  Nystrom  quadrature  for  all  line  integrals  along  dii.  This 
implies  that  one  has  to  interpolate  the  FEM  solution  at  the  Nystrom  points.  To  avoid 
this  extra  interpolation  problem,  we  generate  the  FEM  mesh  such  that  all  the  mesh 
vertices  on  dSl  coincide  with  the  Nystrom  points.  We  therefore  simply  read  off  the  FEM 
nodal  solutions  for  the  Nystrom  quadrature. 

Since  the  coupled  finite  element  and  boundary  integral  approach — together  with 
its  discretization — for  the  adjoint,  incremental  forward,  incremental  adjoint  problems 
is  similar,  we  will  present  only  the  sub-problems  and  the  continuity  condition  at  the 
interface  in  the  next  three  sub-sections. 


7.2.  Adjoint  scattering  problem 


Similar  to  the  forward  scattering  problem,  the  interior  sub-problem  for  the  adjoint 
equation  reads 


VVn  +  k2,num  =  0, 
duin 

+  ikum  =  0, 


dn 


in  hi, 
on  dkl, 


while  the  exterior  problem  reads 

V  Vx  +  k2nuex  =  —K  (U  -  Uobs)  , 


due 


dn 


+  ikuex  =  0, 


in  Rd  \  ft, 
on  d ft, 


lim  r^2 


duey 

dr 


+  ikuex  =  0 


The  interior  and  exterior  solutions  are  matched  by  satisfying  the  following  continuity 
condition  at  the  interface  <9ft: 

um  =  uex,  on  <9ft, 

7.3.  Incremental  forward  scattering  problem 

For  the  incremental  forward  problem,  we  choose  the  interior  sub-problem  as 


V2Um  +  k2nUm  =  - hk 2  (U  +  Uic )  , 
dUm 


dn 


+  ikUm  =  0, 


and  the  exterior  sub-problem  as 

W2Uex  +  k2nUex  =  0, 

+  ikUex  =  0, 

=o, 


dUex 

dn 

dUex 


in  hi, 
on  dSl, 

in  Rd  \  n, 
on  dSl, 


lim  r(d  1)/2  , 

r— >oo  \  Or 
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The  interior  and  exterior  solutions  are  matched  by  satisfying  the  following  continuity 
condition  at  the  interface  <912: 


U in  =  Uex,  on  <912, 


7.f.  Incremental  adjoint  scattering  problem 

Similar  the  the  adjoint  scattering  problem,  the  interior  sub-problem  reads 


V2uin  +  k2nuin 

duin  ,in 
+  iku 


<9n 


while  the  exterior  problem  reads 


— hk2u , 

tp, 


in  12, 
on  <912, 


V2uex  +  k2nuex  =  -KU,  in  Rd  \  12, 

duex 

— - b  ikuex  =  ip,  on  Oil, 

an 

lirn  r{d~1)/2  ( ^  =  0. 

r— >00  \  or  I 

The  interior  and  exterior  solutions  are  matched  by  satisfying  the  following  continuity 
condition  at  the  interface  <912: 


u  =  u 


on  dSl, 


For  numerical  results,  a  second  order  FEM  on  an  unstructured  triangular  mesh  is 
used  for  the  interior  sub-problem,  and  a  Nystrom  method  with  240  equally  distributed 
(in  t)  points  is  used  for  the  exterior  sub-problem.  The  wave  number  is  chosen  to  be 
k  =  10,  while  the  number  of  trigonometric  polynomials  is  60  by  taking  M  =  30.  For 
brevity,  only  results  for  pointwise  observation  are  presented  since  those  for  continuous 
observation  are  similar.  The  observational  data  Uobs  is  synthesized  at  Nohs  =  31  points 
equally  distributed  in  the  interval  y  G  [—10,10]  and  at  x  =  b  =  —10,  unless  otherwise 
stated.  For  the  sake  of  convenience,  the  following  simple  inhomogeneity  [23]  is  used 


n(r) 


1  +  c(l  -  r4)2  r  <  1 

1  r  >  1 


(27) 


where  c  is  some  scalar  constant,  in  particular,  we  choose  c  =  0.5  for  the  synthesization. 
As  a  result,  we  can  choose  12  as  the  unit  circle.  Finally,  the  incident  wave  is  assumed  to 
be  of  the  form  Ulc  =  elkx . 

Our  goal  is  to  numerically  show  that  the  Hessian  is  compact  for  any  bounded 
inhomogeneity  n.  However,  for  convenience,  we  choose  n  of  the  form  (27),  and  in 
particular,  we  choose  to  study  the  discrete  Hessian  at  various  values  of  c.  Numerically, 
we  are  able  to  examine  the  necessary  condition  for  the  Hessian  operator  to  be  compact 
(and  hence  the  ill-posedness  of  the  inverse  problem),  namely,  the  convergence  to  zero 


Hessian  Analysis  for  Inverse  Medium  Acoustic  Scattering 


15 


of  the  Hessian  eigenvalues.  However,  even  in  this  case,  it  is  impossible  to  study  all  the 
eigenvalues  since  they  are  countably  infinite.  We  will  therefore  resort  to  investigating  a 
small  dominant  part  of  the  spectrum,  from  which  we  draw  conclusions.  To  the  rest  of  the 
this  section,  we  “measure”  the  degree  of  ill-posedness  by  the  magnitude  of  eigenvalues. 
For  example,  given  two  ill-posed  inverse  problems,  i.e.,  the  Hessian  eigenvalues  decay  to 
zero,  we  say  one  problem  is  more  ill-posed  than  another  if  the  eigenvalues  of  the  former 
are  smaller  than  those  of  the  latter  at  the  same  indices. 

A  second  order  triangular  mesh  with  2738  elements  and  5597  nodes  is  generated, 
which  permits  us  to  represent  the  refractive  index  as 

5597 

=  ^  ^  rimfm-i 

m=l 

where  are  the  nodal  finite  element  basis  functions.  The  continuous  optimization 
variable  n  has  been  cast  into  5597  discrete  nodal  unknowns  nm ,  and  hence  the  Hessian 
is  an  5597  x  5597  matrix.  The  real  and  imaginary  parts  of  the  forward  and  adjoint 
solutions  at  c  =  1  are  shown  in  Figure  1. 

Away  from  the  optimal  solution,  i.e.,  at  c  =  0.5,  the  Hessian  may  not  be 
(semi-)  positive  definite,  and  for  this  reason,  we  will  present  only  the  eigenvalue 
magnitudes.  Figure  2  shows  the  first  1000  eigenvalues  that  are  largest  in  magnitudes 
for  c  =  {1,0.4,0.499,0.5}.  As  can  be  seen,  the  eigenvalues  decay  exponentially  at  the 
optimal  inhomogeneity,  but  the  decay  rate  is  rather  slow  otherwise.  Moreover,  closer 
to  the  optimal  inhomogeneity,  the  eigenvalue  is  smaller  for  a  same  index,  indicating  the 
increasing  ill-posedness  of  the  inverse  problem  as  the  optimal  solution  is  approached. 
Note  that  at  the  optimal  solution,  the  full  Hessian  collapses  to  the  Gauss-Newton  part, 
i.e.,  Hi,  since  "H2  =  0. 

Next,  we  keep  c  =  1  fixed,  but  allow  the  wave  number  k  to  change.  Figure  3 
shows  the  first  1000  eigenvalues  that  are  largest  in  magnitudes  for  k  =  {10,5,1,0.1}. 
It  can  be  observed,  as  the  wave  number  decreases,  so  do  the  Hessian  eigenvalues.  This 
is  expected  since  intuitively  the  larger  the  wave  number  is,  the  easier  the  detection  of 
the  inhomogeneity,  and  hence  the  problem  is  less  ill-posed.  As  can  also  be  seen,  the 
asymptotic  decay  rate  seems  to  be  similar  for  all  cases. 

We  now  keep  c  =  1  and  k  —  10  fixed,  but  let  the  observation  radius  b  vary. 
We  present  in  Figure  4  the  first  1000  eigenvalues  that  are  largest  in  magnitudes 
for  b  =  {1, 10, 100, 10000}.  As  the  observations  are  taken  further  away  from  the 
inhomogeneity  region,  the  Hessian  eigenvalues  are  smaller.  Again,  for  all  cases,  the 
eigenvalues  decay  to  zero,  indicating  the  compactness  of  the  Hessian  operator.  The  result 
suggests  that  observations  should  be  carried  out  as  close  as  possible  to  the  inhomogeneity 
for  the  inverse  problem  to  be  less  ill-posed. 

In  order  to  study  the  affect  of  observations  on  the  ill-posedness  of  the  inverse 
problem,  we  fix  c  =  1,  k  =  10,  and  b  =  —10,  and  let  Nobs  points  be  equally  distributed 
in  the  interval  y  G  [—100, 100].  Figure  5  shows  the  first  1000  eigenvalues  that  are  largest 
in  magnitudes  for  Nobs  =  {1,51, 101, 1001}.  One  can  observe  that  as  more  observation 
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Figure  1.  Real  and  imaginary  parts  of  the  forward  and  adjoint  solutions  for  n  given 
in  (27)  with  c  =  1. 


Figure  2.  Magnitudes  of  the  first  1000  eigenvalues  of  the  Hessian  at 
n  (c  =  {1,  0.4, 0.499,  0.5})  given  in  (27). 
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Figure  3.  Magnitudes  of  the  first  1000  eigenvalues  of  the  Hessian  for  k  =  {10,  5, 1, 0.1}. 


Figure  4.  Magnitudes  of  the  first  1000  eigenvalues  of  the  Hessian  for  b  = 
{1,10,100,10000}. 


points  are  added,  the  inverse  problem  is  less  ill-posed  since  the  eigenvalues  increase. 
That  is,  as  more  information  about  the  inhomogeneity  is  available,  the  problem  of 
reconstructing  it  is  more  well-posed,  agreeing  with  our  intuition. 

Finally,  we  study  the  dependence  on  mesh  refinement  of  the  Gauss-Newton  Hessian 
dominant  spectrum.  Figure  6  shows  the  first  100  dominant  eigenvalues  for  three  different 
mesh  sizes  h  =  {0.1, 0.05,  0.025}.  As  can  be  observed,  the  dominant  part  of  the  spectrum 
is  numerically  independent  of  the  mesh  size.  This  result  is  consistent  with  the  numerical 
results  in  the  first  part  of  our  work  [19]  in  which  we  numerically  show  that  the  Gauss- 
Newton  Hessian  dominant  spectrum  is  independent  of  the  mesh  size  regardless  of  the 
shape. 
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Figure  5.  Magnitudes  of  the  first  1000  eigenvalues  of  the  Hessian  for  Nobs  = 
{1,51,101,1001}. 


Figure  6.  The  first  100  dominant  eigenvalues  of  the  Hessian  for  h  =  {0.1,  0.05, 0.025} 
with  c  =  0.5. 


8.  Conclusions 

We  have  analyzed  the  Hessian  stemming  from  the  inverse  problem  of  scattering  of 
acoustic  waves  due  to  bounded  inhomogeneity.  Unlike  our  companion  paper  on  inverse 
shape  scattering  problems  [19]  for  which  only  the  Gauss-Newton  Hessian  is  compact, 
the  full  Hessian  operator  has  been  shown  to  be  compact  for  inverse  medium  scattering 
problems.  Our  analysis  starts  with  a  study  on  the  smoothness  of  the  scattering 
solution  based  on  the  Newton  potential  theory  and  the  Riesz-Fredholm  framework. 
Then,  together  with  compact  embeddings  in  Holder  and  Sobolev  spaces,  we  are  able 
to  prove  the  compactness  of  the  Hessian  operator  in  both  Holder  space  and  Sobolev 
space  settings,  and  for  both  two  and  three  dimensions.  Our  theoretical  results  have 
been  validated  numerically  in  several  scenarios. 
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